clear
cd "${PATH}/data"
**************************************************************************************************

*************************************************
// DEFINE LOCALS
*************************************************

// SCC VALUES
local list_SCC = "31.71 50 100"


**************************************************************************************************

*************************************************
// RUN REALLOCATION ALGORITHM (FIGURES 4 AND 5)
*************************************************
*NOTE: As loops are time intensive, it is recommended to run one value for the SCC or alternatively one of the reallocation algorithms at a time


** LOOP OVER DIFFERNET VALUES OF SCC
	
forvalues n=1/3 { 
	local SCC_ele: word `n' of `list_SCC'
	
	global SCC = `SCC_ele'
	global SCC_n = `n'
	
	
	** MAIN REALLOCATION CODE 
	do "${PATH}/programs/2_PVoutput_analysis3_10KW_V3.do" // FIGURE 5

	** ROBUSTNESS1: SOLAR RANKING 
	do "${PATH}/programs/2_PVoutput_analysis3_10KW_V3_robust1.do"
	
	** ROBUSTNESS2: GERMANY WIDE DISPATCH
	do "${PATH}/programs/2_PVoutput_analysis3_10KW_v3_robust2.do"

}


*************************************************
// COMBINE GRAPHS FOR PAPER: 
*************************************************
// SCC VALUE
local list_SCC = "31.71 50 100"

// MAIN PLOT WITH RATIOS (FIGURE 4)

forvalues n=1/3 {
local SCC_ele: word `n' of `list_SCC'

use  "PVoutput_analysis3_10KW_v3_SCC`n'.dta", clear

gen ratios_orig = 100 * (v_alloc / v_bench - 1)

gen ratios_sd_plus_1_plus_1 = 100 * (v_alloc_sd_p_1_p_1 / v_bench_sd_p_1_p_1 - 1)
gen ratios_sd_minus_1_minus_1 = 100 * (v_alloc_sd_m_1_m_1 / v_bench_sd_m_1_m_1 - 1)
gen ratios_sd_minus_1_plus_1 = 100 * (v_alloc_sd_m_1_p_1 / v_bench_sd_m_1_p_1 - 1)
gen ratios_sd_plus_1_minus_1 = 100 * (v_alloc_sd_p_1_m_1 / v_bench_sd_p_1_m_1 - 1)
// CHECK WHICH OF THE ERROR BANDS GENERATES LARGEST AND LOWEST VALUE:
egen ratios_sd_min = rowmin(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)
egen ratios_sd_max = rowmax(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)

keep gamma ratios_orig ratios_sd_min ratios_sd_max


// MERGE WITH ROBUST 1: ALLOCATION ONLY BASED ON SUNSHINE RADIATION
merge 1:1 gamma using "PVoutput_analysis3_10KW_v3_SCC`n'_robust1.dta"
gen ratios_robust1 = 100 * (v_alloc / v_bench - 1)

gen ratios_sd_plus_1_plus_1 = 100 * (v_alloc_sd_p_1_p_1 / v_bench_sd_p_1_p_1 - 1)
gen ratios_sd_minus_1_minus_1 = 100 * (v_alloc_sd_m_1_m_1 / v_bench_sd_m_1_m_1 - 1)
gen ratios_sd_minus_1_plus_1 = 100 * (v_alloc_sd_m_1_p_1 / v_bench_sd_m_1_p_1 - 1)
gen ratios_sd_plus_1_minus_1 = 100 * (v_alloc_sd_p_1_m_1 / v_bench_sd_p_1_m_1 - 1)
// CHECK WHICH OF THE ERROR BANDS GENERATES LARGEST AND LOWEST VALUE:
egen ratios_sd_min_robust1 = rowmin(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)
egen ratios_sd_max_robust1 = rowmax(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)

drop ratios_sd_plus* ratios_sd_minus*
keep gamma ratios*



// MERGE WITH ROBUST 2: ALLOCATION with common merit-order
merge 1:1 gamma using "PVoutput_analysis3_10KW_v3_SCC`n'_robust2.dta"
gen ratios_robust2 = 100 * (v_alloc / v_bench - 1)

gen ratios_sd_plus_1_plus_1 = 100 * (v_alloc_sd_p_1_p_1 / v_bench_sd_p_1_p_1 - 1)
gen ratios_sd_minus_1_minus_1 = 100 * (v_alloc_sd_m_1_m_1 / v_bench_sd_m_1_m_1 - 1)
gen ratios_sd_minus_1_plus_1 = 100 * (v_alloc_sd_m_1_p_1 / v_bench_sd_m_1_p_1 - 1)
gen ratios_sd_plus_1_minus_1 = 100 * (v_alloc_sd_p_1_m_1 / v_bench_sd_p_1_m_1 - 1)
// CHECK WHICH OF THE ERROR BANDS GENERATES LARGEST AND LOWEST VALUE:
egen ratios_sd_min_robust2 = rowmin(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)
egen ratios_sd_max_robust2 = rowmax(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)

drop ratios_sd_plus* ratios_sd_minus*
keep gamma ratios*


tw 	line ratios_orig gamma, lw(medthick) col(navy) || ///
	rarea ratios_sd_max ratios_sd_min gamma, color(gray%20) || ///
	line ratios_robust1 gamma, lw(medthick) lpattern(".-") col(dkgreen) || ///
	rarea ratios_sd_max_robust1 ratios_sd_min_robust1 gamma, color(gray%20) || ///
	line ratios_robust2 gamma, lw(medthick) lpattern("__") col(maroon) || ///
	rarea ratios_sd_max_robust2 ratios_sd_min_robust2 gamma, color(gray%20) ///
	legend(label(1 "TSOs ranking") label(3 "Irradiation ranking") label(5 "Common merit-order") order(1 3 5)) ///
	ylabel(-16(8)25 , grid) yline(0, lc(gray%80)) ytitle("% of gains") xtitle(`"{fontface "Times":{&gamma}}"') ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(ratios_`n', replace)

graph export "../figures/ratios_combined_SCC`n'.pdf", replace 
}
*************************************************


*************************************************
// VALUES vs. GAINS
*APPEND ALL DATA TOGETHER TO ANALYZE ALLOCATION * 

// MAIN PLOT WITH RATIOS: loop over different values for SCC
forvalues n=1/3 {

if (`n'== 1) {
use  "PVoutput_analysis3_10KW_v3_SCC`n'.dta", clear
gen robust =0
}

else {
append using "PVoutput_analysis3_10KW_v3_SCC`n'.dta", gen(appended)
replace robust = 0 if appended ==1
drop appended
}

append using "PVoutput_analysis3_10KW_v3_SCC`n'_robust1.dta", gen(appended)
replace robust = 1 if appended ==1
drop appended

append using "PVoutput_analysis3_10KW_v3_SCC`n'_robust2.dta" , gen(appended)
replace robust = 2 if appended ==1
drop appended
}


******
sort  gamma robust SCC 
br SCC gamma robust v_bench v_alloc

gen diff_alloc = v_alloc - v_bench
tostring SCC, gen(SCC_str) force u // need as string otherwise indicator

// TOTAL VALUE
tw 	(line v_alloc gamma if robust == 0 & SCC_str == "31.71", lc(navy)) /// 
	(line v_alloc gamma if robust == 1 & SCC_str =="31.71", lp(_)  lc(dkgreen)) ///
	(line v_alloc gamma if robust == 2 & SCC_str =="31.71", lp(.-) lc(maroon)) ///
	(line v_alloc gamma if robust == 0 & SCC_str =="50", lc(navy)) /// 
	(line v_alloc gamma if robust == 1 & SCC_str =="50", lp(_) lc(dkgreen))  ///
	(line v_alloc gamma if robust == 2 & SCC_str =="50" , lp(.-) lc(maroon)) ///
	(line v_alloc gamma if robust == 0 & SCC_str == "100", lc(navy)) /// 
	(line v_alloc gamma if robust == 1 & SCC_str == "100", lp(_) lc(dkgreen)) ///
	(line v_alloc gamma if robust == 2 & SCC_str == "100" , lp(.-) lc(maroon)), ///
	legend(label(1 "TSOs ranking") label(2 "Irradiation ranking") label(3 "Common merit-order") order(1 2 3) rows(1) symysize(*.8)) ///
	ylabel(1500(1000)4000 , grid)  ///
	text(3800 .37 "SCC = 100 {c 0128}/tCO{subscript:2}",place(e) size(*.7)) ///
	text(2200 .37 "SCC = 50 {c 0128}/tCO{subscript:2}",place(e) size(*.7)) ///
	text(1700 .37 "SCC = 31.71 {c 0128}/tCO{subscript:2}",place(e) size(*.7)) ///
	ytitle("Total value of reallocation [m{c 0128}]") xtitle(`"{fontface "Times":{&gamma}}"') ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(value_realloc, replace)

graph export "../figures/value_reallocation.pdf", replace



**************************************************************************************************
**************************************************************************************************
// ROBUSTNESS: RERUN MAIN REALLOCATION INCLUDING CO-POLLUTANTS (ONLINE APPENDIX FIGURE D.4)

*************************************************
// RUN REALLOCATION ALGORITHM
*************************************************// 

	** MAIN REALLOCATION CODE 
	do "${PATH}/programs/2_PVoutput_analysis3_10KW_v4_lpollution.do" // Figure D.4, Panels b-d
	
	** ROBUSTNESS1: SOLAR RANKING 
	do "${PATH}/programs/2_PVoutput_analysis3_10KW_v4_lpollution_robust1.do"
	
	** ROBUSTNESS2: GERMANY WIDE DISPATCH
	do "${PATH}/programs/2_PVoutput_analysis3_10KW_v4_lpollution_robust2.do"



*************************************************
// COMBINE GRAPHS FOR PAPER: Figure D.4, Panel a
*************************************************
// ORIGINAL ALLOCATION
use "PVoutput_analysis3_10KW_v4_lpollution_SCC1", clear
gen ratios_orig = 100 * (v_alloc / v_bench - 1)
keep gamma ratios_orig

// MERGE WITH ROBUST 1: ALLOCATION ONLY BASED ON SUNSHINE RADIATION
merge 1:1 gamma using "PVoutput_analysis3_10KW_v4_lpollution_SCC1_robust1.dta"

gen ratios_robust1 = 100 * (v_alloc / v_bench - 1)
keep gamma ratios*


// MERGE WITH ROBUST 2: ALLOCATION with common merit-order
merge 1:1 gamma using "PVoutput_analysis3_10KW_v4_lpollution_SCC1_robust2.dta"

gen ratios_robust2 = 100 * (v_alloc / v_bench - 1)
keep gamma ratios*


tw 	line ratios_orig gamma, lw(medthick) col(navy) || ///
	line ratios_robust1 gamma, lw(medthick) lpattern(".-") col(dkgreen) || ///
	line ratios_robust2 gamma, lw(medthick) lpattern("__") col(maroon) ///
		legend(label(1 "TSOs ranking") label(2 "Irradiation ranking") label(3 "Common merit-order")) ylabel(-10(10)30 , grid) ytitle("% of gains") xtitle(`"{fontface "Times":{&gamma}}"') ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(ratios, replace)

graph export "../figures/ratios_v4_lpollution_SCC1_combined.pdf", replace 



